MakeIndex <- function(q1, q0) {
	count.stats <- 1:(q1  + q0)
	index.treat <- combn(count.stats, q1)
	index.cont <- sapply(1:ncol(index.treat), function(j) count.stats[! count.stats %in% index.treat[, j]])	
	index <- cbind(t(index.treat), t(index.cont))
	index
}

OrdTest <- function(stats, k, index, q1, q0) {
	q <- q1 + q0
	M <- nrow(index)
	r.stats <- matrix(stats[index[1:M, ]], M, q)
	perm.stats <- rowMeans(r.stats[, 1:q1]) - rowMeans(r.stats[, (q1+1):q])
	orig.stat <- perm.stats[1]
	perm.stats <- sort(perm.stats)
	as.numeric(orig.stat > perm.stats[k])
}

IMTest <- function(stats, alpha = .05, q1, q0) {
	q <- q1 + q0
	t.stat <- (mean(stats[1:q1]) - mean(stats[(q1+1):q]))/sqrt(var(stats[1:q1])/q1 + var(stats[(q1+1):q])/q0)
	as.numeric(t.stat > qt(1-alpha, min(q1, q0) - 1))
}

d <- read.csv("bloometal13.csv")

index <- MakeIndex(11, 6)

# .008 is value from Table 1 for 5% test
k <- ceiling((1-.008)*dim(index)[1])

# Reject null at 5%
OrdTest(d[, 2], k, index, 11, 6)